In [1]:
# necessary packages #

#using Pkg
#Pkg.add("Distances")
using Distributions
using Random
using Distances
using LinearAlgebra
using SparseArrays
using IterativeSolvers
using ProgressMeter
In [2]:
include("util.j")
Out[2]:
getAD (generic function with 1 method)
In [24]:
# unnecessary packages #

#using Pkg
#Pkg.add("UnicodePlots")
using UnicodePlots   # check the structure of the sparse matrix
using BenchmarkTools

using StatsPlots
using MCMCChains
using PrettyTables
In [4]:
#using Pkg
#Pkg.add("ProgressMeter");
In [5]:
# Set the parameters for SLMC model #

N = 1200 # sample size
N1 = 1000; N2 = 1000;
q = 2; p = 2; K = 2
Σ = [0.3 0.1
     0.1 0.2];
β = [1.0 -1.0
     -5.0 2.0];
ϕ1 = 3.0; ϕ2 = 30.0; ν1 = 0.5; ν2 = 0.5; # parameter for the independent F
Λ = [3.0 -4.0
     -2.0 3.0];

# priors #
μβ = fill(0.0, p, q); Vβ =[[100.0 0.0]; [0.0 100.0]];
μΛ = fill(0.0, K, q); VΛ =[[100.0 0.0]; [0.0 100.0]];
νΣ = 2 * q; ΨΣ = [[1.0 0.0]; [0.0 1.0]];
In [6]:
# Generate simulation data #

Random.seed!(1234);
coords = rand(2, N);                                          # random location over unit square
X = hcat(fill(1, (N,)), rand(N));                             # design matrix
D = pairwise(Euclidean(), coords, dims = 2);                  # distance matrix
ρ1 = exp.(-ϕ1 * D); ρ2 = exp.(-ϕ2 * D);                       # covariance matrix
ω = [rand(MvNormal(ρ1), 1) rand(MvNormal(ρ2), 1)] * Λ; # latent process
Y = X * β + ω + transpose(rand(MvNormal(Σ), N));              # response matrix
In [7]:
# Some data preparations #

ordx = sortperm(coords[1, :]);                                # sort order based on the first coordinates
X_ord = X[ordx, :]; Y_ord = Y[ordx, :]; ω_ord = ω[ordx, :];   # sorted data
ω_incp_obs = ω_ord + fill(1.0, (N, 1)) * transpose(β[1, :]);  # latent process + intercept
coords_ord = coords[:, ordx];
S1_ind = sample(1:N, N1, replace = false, ordered = true);    # observed location index for 1st response
S2_ind = sample(1:N, N2, replace = false, ordered = true);    # observed location index for 2nd response
S = sort(union(S1_ind, S2_ind));                              # observed index set
M1_ind = setdiff(S, S1_ind);                                  # in S not in S1
M2_ind = setdiff(S, S2_ind);                                  # in S not in S2 
obs_ind = vcat(S1_ind, S2_ind .+ N);              # index of the observed location for all response among N locations
perm_ind = sortperm(vcat(S1_ind, S2_ind));                    # the vector of the permutation 

v1 = zeros(N); v1[S1_ind] .= 1;
v2 = zeros(N); v2[S2_ind] .= 1;
index_S = (2^0 * v1 + 2^1 * v2);                              # build index indicating which response are observed
M1_Sind = findall(x -> x == 2^1, index_S[S]);                 # index of M1 among S
M2_Sind = findall(x -> x == 2^0, index_S[S]);                 # index of M2 among S

m = 10; n = length(S);                                        # number of nearest neighbor                       
NN = BuildNN(coords_ord[:, S], m, 1.0);                            # build nearest neighbor 
nnIndx_col = vcat(NN.nnIndx, 1:n);                            # the index of columns
nnIndx_row = zeros(Int64, 0);                                               
for i in 2:m
    nnIndx_row = vcat(nnIndx_row, fill(i, i-1))
end
nnIndx_row = vcat(nnIndx_row, repeat((m + 1):n, inner = m), 1:n);  # the index of rows

Lβ = cholesky(Vβ); LΛ = cholesky(VΛ);
In [8]:
# check the plot of the data 
using RCall
@rput ω_ord
@rput coords_ord
@rput S
R"""
library(MBA)
library(classInt)
library(RColorBrewer)
library(sp)
library(coda)
library(spBayes)
library(fields)

h <- 12
surf.raw1 <- mba.surf(cbind(t(coords_ord[, S]), ω_ord[S, 1]), no.X = 300, no.Y = 300, 
              exten = TRUE, sp = TRUE, h = h)$xyz.est
surf.raw2 <- mba.surf(cbind(t(coords_ord[, S]), ω_ord[S, 2]), no.X = 300, no.Y = 300, 
              exten = TRUE, sp = TRUE, h = h)$xyz.est
surf.brks <- classIntervals(c(surf.raw1[["z"]], surf.raw2[["z"]]), 500, 'pretty')$brks
col.pal <- colorRampPalette(brewer.pal(11,'RdBu')[11:1])
xlim <- c(0, 1.13)

zlim <- range(c(surf.raw1[["z"]], surf.raw2[["z"]]))

# size for the mapping of w               
width <- 360
height <- 360
pointsize <- 16

png(paste("./pics/map-w1-true-Keqq.png", sep = ""), 
    width = width, height = height, pointsize = pointsize, family = "Courier")
par(mfrow = c(1, 1))
##Obs
i <- as.image.SpatialGridDataFrame(surf.raw1)
plot(t(coords_ord[, S]), typ="n", cex=0.5, xlim=xlim, axes=FALSE, ylab="y", xlab="x") 
     #main = "true")
axis(2, las=1)
axis(1)
image.plot(i, add=TRUE, col=rev(col.pal(length(surf.brks)-1)), zlim=zlim)
dev.off()

png(paste("./pics/map-w2-true-Keqq.png", sep = ""), 
    width = width, height = height, pointsize = pointsize, family = "Courier")
par(mfrow = c(1, 1))
##Obs
i <- as.image.SpatialGridDataFrame(surf.raw2)
plot(t(coords_ord[, S]), typ="n", cex=0.5, xlim=xlim, axes=FALSE, ylab="y", xlab="x") 
     #main = "true")
axis(2, las=1)
axis(1)
image.plot(i, add=TRUE, col=rev(col.pal(length(surf.brks)-1)), zlim=zlim)
dev.off()

"""

┌ Warning: Error requiring AxisArrays from RCall:
│ LoadError: ArgumentError: Package RCall does not have AxisArrays in its dependencies:
│ - If you have RCall checked out for development and have
│   added AxisArrays as a dependency but haven't updated your primary
│   environment's manifest file, try `Pkg.resolve()`.
│ - Otherwise you may need to report an issue with RCall
│ Stacktrace:
│  [1] require(::Module, ::Symbol) at ./loading.jl:836
│  [2] include at ./boot.jl:326 [inlined]
│  [3] include_relative(::Module, ::String) at ./loading.jl:1038
│  [4] include at ./sysimg.jl:29 [inlined]
│  [5] include(::String) at /home/lu/.julia/packages/RCall/ffM0W/src/RCall.jl:2
│  [6] top-level scope at /home/lu/.julia/packages/RCall/ffM0W/src/setup.jl:179
│  [7] eval at ./boot.jl:328 [inlined]
│  [8] eval at /home/lu/.julia/packages/RCall/ffM0W/src/RCall.jl:2 [inlined]
│  [9] (::getfield(RCall, Symbol("##5#8")))() at /home/lu/.julia/packages/Requires/9Jse8/src/require.jl:67
│  [10] err(::getfield(RCall, Symbol("##5#8")), ::Module, ::String) at /home/lu/.julia/packages/Requires/9Jse8/src/require.jl:38
│  [11] #4 at /home/lu/.julia/packages/Requires/9Jse8/src/require.jl:66 [inlined]
│  [12] withpath(::getfield(RCall, Symbol("##4#7")), ::String) at /home/lu/.julia/packages/Requires/9Jse8/src/require.jl:28
│  [13] #3 at /home/lu/.julia/packages/Requires/9Jse8/src/require.jl:65 [inlined]
│  [14] listenpkg(::getfield(RCall, Symbol("##3#6")), ::Base.PkgId) at /home/lu/.julia/packages/Requires/9Jse8/src/require.jl:13
│  [15] macro expansion at /home/lu/.julia/packages/Requires/9Jse8/src/require.jl:64 [inlined]
│  [16] __init__() at /home/lu/.julia/packages/RCall/ffM0W/src/setup.jl:178
│  [17] _include_from_serialized(::String, ::Array{Any,1}) at ./loading.jl:633
│  [18] _require_search_from_serialized(::Base.PkgId, ::String) at ./loading.jl:713
│  [19] _require(::Base.PkgId) at ./loading.jl:937
│  [20] require(::Base.PkgId) at ./loading.jl:858
│  [21] require(::Module, ::Symbol) at ./loading.jl:853
│  [22] top-level scope at In[8]:1
│  [23] eval at ./boot.jl:328 [inlined]
│  [24] softscope_include_string(::Module, ::String, ::String) at /home/lu/.julia/packages/SoftGlobalScope/dLfaZ/src/SoftGlobalScope.jl:218
│  [25] execute_request(::ZMQ.Socket, ::IJulia.Msg) at /home/lu/.julia/packages/IJulia/4UizY/src/execute_request.jl:67
│  [26] #invokelatest#1 at ./essentials.jl:742 [inlined]
│  [27] invokelatest at ./essentials.jl:741 [inlined]
│  [28] eventloop(::ZMQ.Socket) at /home/lu/.julia/packages/IJulia/4UizY/src/eventloop.jl:8
│  [29] (::getfield(IJulia, Symbol("##15#18")))() at ./task.jl:259
│ in expression starting at /home/lu/.julia/packages/RCall/ffM0W/src/convert/axisarray.jl:1
└ @ Requires /home/lu/.julia/packages/Requires/9Jse8/src/require.jl:40
┌ Warning: RCall.jl: Loading required package: magic
│ Loading required package: abind
│ Loading required package: Formula
│ Loading required package: Matrix
└ @ RCall /home/lu/.julia/packages/RCall/ffM0W/src/io.jl:113
┌ Warning: RCall.jl: Loading required package: spam
│ Loading required package: dotCall64
│ Loading required package: grid
│ Spam version 2.2-2 (2019-03-07) is loaded.
│ Type 'help( Spam)' or 'demo( spam)' for a short introduction 
│ and overview of this package.
│ Help for individual functions is also obtained by adding the
│ suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
│ 
│ Attaching package: ‘spam’
│ 
│ The following object is masked from ‘package:Matrix’:
│ 
│     det
│ 
│ The following objects are masked from ‘package:base’:
│ 
│     backsolve, forwardsolve
│ 
│ Loading required package: maps
│ See www.image.ucar.edu/~nychka/Fields for
│  a vignette and other supplements. 
└ @ RCall /home/lu/.julia/packages/RCall/ffM0W/src/io.jl:113
Out[8]:
RObject{IntSxp}
null device 
          1
In [9]:
# check neighbors
index = 660
gr()
scatter(coords_ord[1, S], coords_ord[2, S], xlim = (0, 1), ylim = (0, 1), color = :white)
scatter!(coords_ord[1, S[1:index]], coords_ord[2, S[1:index]], color = :grey)
scatter!(coords_ord[1, S[index:index]], coords_ord[2, S[index:index]], color = :red)
scatter!(coords_ord[1, S[NN.nnIndx[NN.nnIndxLU[index - 1]:(NN.nnIndxLU[index] - 1)]]], 
    coords_ord[2, S[NN.nnIndx[NN.nnIndxLU[index - 1]:(NN.nnIndxLU[index] - 1)]]], color = :orange)
Out[9]:
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 y1 y2 y3 y4
In [10]:
# preallocation #

#F = Array{Float64,2}(undef, n , 3);                           # preallocate the matrix F
μ_m1 = Array{Float64, 2}(undef, length(M1_ind), q);
μ_m2 = Array{Float64, 2}(undef, length(M2_ind), q);
nIndx = length(NN.nnIndx);
A1 = Array{Float64}(undef, nIndx); D1 = Array{Float64}(undef, n);
A2 = Array{Float64}(undef, nIndx); D2 = Array{Float64}(undef, n);
Ystar = vcat(Y_ord[S, :], Lβ.L \ μβ, LΛ.L \ μΛ);             # will be updated after imputing missing response
Xstar = vcat([X_ord[S, :] fill(0.0, n, K)], [inv(Lβ.L) fill(0.0, p, K)], 
    [fill(0.0, K, p) inv(LΛ.L)]);      
Ψstar = fill(0.0, q, q); νstar = νΣ + n;
μγstar = vcat(μβ, μΛ); Vγstar = fill(0.0, p + K, p + K);
Y_Xm = fill(0.0, n, q); invVγstar = fill(0.0, p + K, p + K);

MCMC sampling algorithm Q1: priors for $\nu_i$ Q2: $\phi_i$ may not be consistant, since the order can change

In [11]:
# Preallocation for MCMC samples and Initalization #
N_sam = 20000;
γ_sam = Array{Float64, 3}(undef, p + K, q, N_sam + 1);
Σ_sam = Array{Float64, 3}(undef, q, q, N_sam + 1);
F_sam = Array{Float64, 3}(undef, n, K, N_sam);
Y_m1_sam = Array{Float64, 2}(undef, length(M1_ind), N_sam);
Y_m2_sam = Array{Float64, 2}(undef, length(M2_ind), N_sam);

ϕ1_0 = ϕ1; ϕ2_0 = ϕ2; 
γ_sam[:, :, 1] = vcat([[0.0 0.0]; [0.0 0.0]], [[1.0 0.0]; [0.0 1.0]]);
Σ_sam[:, :, 1] = [[0.5 0.1]; [0.1 0.5]];

# first consider SLMC with fixed hyperparameter set {ψk} #
# Build the matrix Vk #
getAD(coords_ord[:, S], NN.nnIndx, NN.nnDist, NN.nnIndxLU, ϕ1_0, 0.5, A1, D1);
getAD(coords_ord[:, S], NN.nnIndx, NN.nnDist, NN.nnIndxLU, ϕ2_0, 0.5, A2, D2);
In [12]:
# for loop for MCMC chain #
prog = Progress(N_sam, 1, "Computing initial pass...", 50)
for l in 1:N_sam
    # Build the matrix D_Sigma_o^{1/2} #
    Dic_diag = Dict(2^0 => sparse(1I, 1, 1) * (1 / sqrt(Σ_sam[:, :, l][1, 1])), 
        2^1 => sparse(1I, 1, 1) * (1 / sqrt(Σ_sam[:, :, l][2, 2])), 
        (2^0 + 2^1)=> sparse(sqrt(inv(Σ_sam[:, :, l]))));
    invD = blockdiag([Dic_diag[i] for i in index_S if i > 0]...);
                    
    # first consider SLMC with fixed hyperparameter set {ψk} #
    # Build the matrix Vk #
    #getAD(coords_ord[:, S], NN.nnIndx, NN.nnDist, NN.nnIndxLU, ρ1_0, 0.5, A1, D1);
    #getAD(coords_ord[:, S], NN.nnIndx, NN.nnDist, NN.nnIndxLU, ρ2_0, 0.5, A2, D2);
    #getAD(coords_ord[:, S], NN.nnIndx, NN.nnDist, NN.nnIndxLU, ρ3_0, 0.5, A3, D3);

    # Build Ytilde Xtilde
    Ytilde = vcat(invD * vcat(Y_ord[S1_ind, 1] - X_ord[S1_ind, :] * γ_sam[1:p, 1, l], 
                            Y_ord[S2_ind, 2] - X_ord[S2_ind, :] * γ_sam[1:p, 2, l])[perm_ind], zeros(K * n));
    Xtilde = vcat(
             invD * kron(sparse(transpose(γ_sam[(p + 1):(p + K), :, l])), 
                            sparse(1:N, 1:N, ones(N)))[obs_ind, 
                            vcat(S, S .+ N)][perm_ind, :],
             blockdiag(
             Diagonal(1 ./ sqrt.(D1)) * sparse(nnIndx_row, nnIndx_col, vcat(-A1, ones(n))),
             Diagonal(1 ./ sqrt.(D2)) * sparse(nnIndx_row, nnIndx_col, vcat(-A2, ones(n)))))
                
    # use LSMR to generate sample of F #       
    nsam = length(Ytilde);
    F_sam[:, :, l] = reshape(lsmr(Xtilde, Ytilde + rand(Normal(), nsam)), :, K);
                    
    # impute missing response  over S#
    Xstar[1:n, (p + 1):(p + K)] = F_sam[:, :, l];        # update matrix Xstar with F
    mul!(μ_m1, Xstar[M1_Sind, :], γ_sam[:, :, l]);
    mul!(μ_m2, Xstar[M2_Sind, :], γ_sam[:, :, l]);

    Y_m1_sam[:, l] = μ_m1[:, 1] + (Σ_sam[1, 2, l] / Σ_sam[2, 2, l]) * 
            (Y_ord[M1_ind, 2] - μ_m1[:, 2]) + 
            rand(Normal(0, sqrt(Σ_sam[1, 1, l] - Σ_sam[1, 2, l]^2 / Σ_sam[2, 2, l])), length(M1_ind));
    Y_m2_sam[:, l] = μ_m2[:, 2] + (Σ_sam[2, 1, l] / Σ_sam[1, 1, l]) * 
            (Y_ord[M2_ind, 1] - μ_m2[:, 1]) + 
            rand(Normal(0, sqrt(Σ_sam[2, 2, l] - Σ_sam[2, 1, l]^2 / Σ_sam[1, 1, l])), length(M2_ind)); # improve?...
                    
    # use MNIW to sample γ Σ #
    Ystar[M1_Sind, 1] = Y_m1_sam[:, l];              # update Ystar with imputed response
    Ystar[M2_Sind, 2] = Y_m2_sam[:, l]; 
    invVγstar = cholesky(Xstar'Xstar);
    mul!(μγstar, transpose(Xstar), Ystar); μγstar = invVγstar.U \ (invVγstar.L \ μγstar);
    Y_Xm = BLAS.gemm('N', 'N', -1.0, Xstar, μγstar) + Ystar;
    mul!(Ψstar, transpose(Y_Xm), Y_Xm); BLAS.axpy!(1.0, ΨΣ, Ψstar);

    Σ_sam[:, :, l + 1] = rand(InverseWishart(νstar, Ψstar), 1)[1];    # sample Σ
    γ_sam[:, :, l + 1] = (invVγstar.U \ reshape(rand(Normal(), (p + K) * q), (p + K), q)) * 
                    cholesky(Σ_sam[:, :, l + 1]).U + μγstar;          # sample γ    
    next!(prog) # monitor the progress
end

Computing initial pass...100%|██████████████████████████████████████████████████| Time: 1:11:37:50

MCMC Chain check

In [13]:
β_pos_sam = Array{Float64, 3}(undef, N_sam + 1, p * q, 1);
β_pos_sam[:, :, 1] = hcat(γ_sam[1, 1, :], γ_sam[1, 2, :], γ_sam[2, 1, :], γ_sam[2, 2, :]);
β_chain = Chains(β_pos_sam);
pβ = plot(β_chain)
Out[13]:
0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -0.5 0.0 0.5 1.0 Param1 Iteration Sample value -1.0 -0.5 0.0 0.5 1.0 0.00 0.25 0.50 0.75 1.00 Param1 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -1.0 -0.5 0.0 0.5 1.0 1.5 Param2 Iteration Sample value -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 0.0 0.2 0.4 0.6 0.8 Param2 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -5 -4 -3 -2 -1 0 Param3 Iteration Sample value -5 -4 -3 -2 -1 0 0.0 0.5 1.0 1.5 2.0 2.5 Param3 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0.0 0.5 1.0 1.5 2.0 2.5 Param4 Iteration Sample value 0.0 0.5 1.0 1.5 2.0 2.5 0.0 0.5 1.0 1.5 2.0 Param4 Sample value Density
In [14]:
Λ_pos_sam = Array{Float64, 3}(undef, N_sam + 1, K * q, 1);
Λ_pos_sam[:, :, 1] = hcat(γ_sam[3, 1, :], γ_sam[3, 2, :], γ_sam[4, 1, :], γ_sam[4, 2, :]);
Λ_chain = Chains(Λ_pos_sam);
pΛ = plot(Λ_chain)
Out[14]:
0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -3 -2 -1 0 1 Param1 Iteration Sample value -4 -3 -2 -1 0 1 0.0 0.2 0.4 0.6 0.8 Param1 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 1 2 3 4 5 Param2 Iteration Sample value 0 1 2 3 4 5 0.0 0.1 0.2 0.3 0.4 0.5 0.6 Param2 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 Param3 Iteration Sample value -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0 1 2 3 Param3 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 1.0 1.5 2.0 2.5 3.0 3.5 Param4 Iteration Sample value 1.0 1.5 2.0 2.5 3.0 3.5 0 1 2 3 Param4 Sample value Density
In [15]:
Σ_pos_sam = Array{Float64, 3}(undef, N_sam + 1, q * q, 1);
Σ_pos_sam[:, :, 1] = hcat(Σ_sam[1, 1, :], Σ_sam[1, 2, :], Σ_sam[2, 1, :], Σ_sam[2, 2, :]);
Σ_chain = Chains(Σ_pos_sam);
pΣ = plot(Σ_chain)
Out[15]:
0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0.3 0.6 0.9 1.2 1.5 Param1 Iteration Sample value 0.3 0.6 0.9 1.2 1.5 0 1 2 3 4 5 Param1 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -1.5 -1.0 -0.5 0.0 Param2 Iteration Sample value -1.5 -1.0 -0.5 0.0 0 2 4 6 8 Param2 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -1.5 -1.0 -0.5 0.0 Param3 Iteration Sample value -1.5 -1.0 -0.5 0.0 0 2 4 6 8 Param3 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 1 2 3 Param4 Iteration Sample value 0 1 2 3 0 2 4 6 8 Param4 Sample value Density
In [25]:
ω_incp_obs_pos_sam = Array{Float64, 3}(undef, n, q, N_sam);
lll = fill(1.0, (n, 1));
for i in 1:N_sam
    ω_incp_obs_pos_sam[:, :, i] = F_sam[:, :, i] * γ_sam[(p + 1):(p + K), :, i + 1] + 
                     lll * transpose(γ_sam[1, :, i + 1]);
end
truncindex = 1;#Integer(trunc(N_sam / 2));
ω_incp_pos_sam = Array{Float64, 3}(undef, N_sam - truncindex  + 1, 3, 1);
ω_incp_pos_sam[:, :, 1] = hcat(ω_incp_obs_pos_sam[1, 1, truncindex:N_sam], 
    ω_incp_obs_pos_sam[1, 2, truncindex:N_sam], ω_incp_obs_pos_sam[200, 1, truncindex:N_sam]);
ω_incp_chain = Chains(ω_incp_pos_sam);
pΣ = plot(ω_incp_chain)
Out[25]:
0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -6 -5 -4 -3 -2 Param1 Iteration Sample value -6 -5 -4 -3 -2 0.0 0.2 0.4 0.6 Param1 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 2 4 6 8 Param2 Iteration Sample value 2 4 6 8 10 0.0 0.1 0.2 0.3 0.4 Param2 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 1 2 3 4 5 Param3 Iteration Sample value 0 1 2 3 4 5 0.0 0.2 0.4 0.6 Param3 Sample value Density

Posterior Inference

In [26]:
N_Inf_burn = Integer(trunc(N_sam / 2));
ω_incp_obs_pos_qt = Array{Float64, 3}(undef, n, q, 3);
for j in 1:q
    for i in 1:n
        ω_incp_obs_pos_qt[i, j, :] = quantile(ω_incp_obs_pos_sam[i, j, N_Inf_burn:N_sam], [0.025, 0.5, 0.975])
    end
end
# count the covarage of 95% CI #
count = fill(0.0, 2);
for j in 1:q
    for i in 1:n
        count[j] = count[j] + 
        ((ω_incp_obs_pos_qt[i, j, 1] < ω_incp_obs[S[i], j]) && 
            (ω_incp_obs_pos_qt[i, j, 3] > ω_incp_obs[S[i], j]))
    end
end
count
Out[26]:
2-element Array{Float64,1}:
 1141.0
 1137.0
In [27]:
count ./ n
Out[27]:
2-element Array{Float64,1}:
 0.9743808710503843
 0.9709649871904356
In [29]:
summary_table = Array{Float64, 2}(undef, (p - 1) * q + q * q, 5);
summary_table[1, :] = vcat(β[2, 1], mean(γ_sam[2, 1, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[2, 1, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[2, :] = vcat(β[2, 2], mean(γ_sam[2, 2, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[2, 2, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[3, :] = vcat(Σ[1, 1], mean(Σ_sam[1, 1, N_Inf_burn:(N_sam + 1)]),
    quantile(Σ_sam[1, 1, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[4, :] = vcat(Σ[1, 2], mean(Σ_sam[1, 2, N_Inf_burn:(N_sam + 1)]),
    quantile(Σ_sam[1, 2, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[5, :] = vcat(Σ[2, 1], mean(Σ_sam[2, 1, N_Inf_burn:(N_sam + 1)]),
    quantile(Σ_sam[2, 1, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[6, :] = vcat(Σ[2, 2], mean(Σ_sam[2, 2, N_Inf_burn:(N_sam + 1)]),
    quantile(Σ_sam[2, 2, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table = round.(summary_table; digits = 3 );
rnames = ["β[2, 1]", "β[2, 2]", "Σ[1, 1]", "Σ[1, 2]", "Σ[2, 1]", "Σ[2, 2]"];
summary_table = [rnames summary_table];
pretty_table(summary_table,  ["" "true" "mean" "median" "2.5%" "97.5%"], markdown)

|         | true |   mean | median |   2.5% |  97.5% |
|---------|------|--------|--------|--------|--------|
| β[2, 1] | -5.0 | -4.977 | -4.976 | -5.265 | -4.694 |
| β[2, 2] |  2.0 |  1.932 |  1.934 |  1.548 |  2.336 |
| Σ[1, 1] |  0.3 |  0.356 |  0.351 |   0.23 |  0.503 |
| Σ[1, 2] |  0.1 |  0.053 |  0.057 | -0.046 |  0.131 |
| Σ[2, 1] |  0.1 |  0.053 |  0.057 | -0.046 |  0.131 |
| Σ[2, 2] |  0.2 |  0.149 |  0.144 |  0.068 |  0.254 |

Plot the latent process + intercept

In [17]:
# check the plot of the data 
using RCall
@rput ω_incp_obs
@rput coords_ord
@rput S
@rput ω_incp_obs_pos_qt
R"""
library(MBA)
library(classInt)
library(RColorBrewer)
library(sp)
library(coda)
library(spBayes)
library(fields)

h <- 12
surf.raw1 <- mba.surf(cbind(t(coords_ord[, S]), ω_incp_obs[S, 1]), no.X = 300, no.Y = 300, 
              exten = TRUE, sp = TRUE, h = h)$xyz.est
surf.raw2 <- mba.surf(cbind(t(coords_ord[, S]), ω_incp_obs[S, 2]), no.X = 300, no.Y = 300, 
              exten = TRUE, sp = TRUE, h = h)$xyz.est

surf.pos1 <- mba.surf(cbind(t(coords_ord[, S]), ω_incp_obs_pos_qt[, 1, 2]), no.X = 300, no.Y = 300, 
              exten = TRUE, sp = TRUE, h = h)$xyz.est
surf.pos2 <- mba.surf(cbind(t(coords_ord[, S]), ω_incp_obs_pos_qt[, 2, 2]), no.X = 300, no.Y = 300, 
              exten = TRUE, sp = TRUE, h = h)$xyz.est

surf.brks <- classIntervals(c(surf.raw1[["z"]], surf.raw2[["z"]]), 500, 'pretty')$brks
col.pal <- colorRampPalette(brewer.pal(11,'RdBu')[11:1])
xlim <- c(0, 1.13)

zlim <- range(c(surf.raw1[["z"]], surf.raw2[["z"]], surf.pos1[["z"]], surf.pos2[["z"]]))

# size for the mapping of w               
width <- 360
height <- 360
pointsize <- 16

png(paste("./pics/map-w1_incp-true_Keqq.png", sep = ""), 
    width = width, height = height, pointsize = pointsize, family = "Courier")
par(mfrow = c(1, 1))
##Obs
i <- as.image.SpatialGridDataFrame(surf.raw1)
plot(t(coords_ord[, S]), typ="n", cex=0.5, xlim=xlim, axes=FALSE, ylab="y", xlab="x") 
     #main = "true")
axis(2, las=1)
axis(1)
image.plot(i, add=TRUE, col=rev(col.pal(length(surf.brks)-1)), zlim=zlim)
dev.off()

png(paste("./pics/map-w2_incp-true_Keqq.png", sep = ""), 
    width = width, height = height, pointsize = pointsize, family = "Courier")
par(mfrow = c(1, 1))
##Obs
i <- as.image.SpatialGridDataFrame(surf.raw2)
plot(t(coords_ord[, S]), typ="n", cex=0.5, xlim=xlim, axes=FALSE, ylab="y", xlab="x") 
     #main = "true")
axis(2, las=1)
axis(1)
image.plot(i, add=TRUE, col=rev(col.pal(length(surf.brks)-1)), zlim=zlim)
dev.off()

png(paste("./pics/map-w1_incp-posm_Keqq.png", sep = ""), 
    width = width, height = height, pointsize = pointsize, family = "Courier")
par(mfrow = c(1, 1))
##Obs
i <- as.image.SpatialGridDataFrame(surf.pos1)
plot(t(coords_ord[, S]), typ="n", cex=0.5, xlim=xlim, axes=FALSE, ylab="y", xlab="x") 
     #main = "true")
axis(2, las=1)
axis(1)
image.plot(i, add=TRUE, col=rev(col.pal(length(surf.brks)-1)), zlim=zlim)
dev.off()

png(paste("./pics/map-w2_incp-posm_Keqq.png", sep = ""), 
    width = width, height = height, pointsize = pointsize, family = "Courier")
par(mfrow = c(1, 1))
##Obs
i <- as.image.SpatialGridDataFrame(surf.pos2)
plot(t(coords_ord[, S]), typ="n", cex=0.5, xlim=xlim, axes=FALSE, ylab="y", xlab="x") 
     #main = "true")
axis(2, las=1)
axis(1)
image.plot(i, add=TRUE, col=rev(col.pal(length(surf.brks)-1)), zlim=zlim)
dev.off()

"""
Out[17]:
RObject{IntSxp}
null device 
          1

![w1] ![w1_pos] ![w2] ![w2_pos]